!
!  FLATBAT     FLAT     Flat field correction utility.
!  Copyright (C) 1997-x  Filip Hroch, Masaryk University, Brno, CZ
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

program flatbat

  use input
  use robustmean

  implicit none

  ! No. of image dimensions
  integer, parameter :: DIM = 2

  ! Version:
  Character(len=*),parameter :: version = 'FLATBAT Version 2.0, (C)1997-2010 F.Hroch (hroch@physics.muni.cz)'

  ! The flatbat assumes that output image has in intensity
  ! within range minvalue < x < maxvalue.
  ! The minvalue always is 0 for measured intensities (?)
  ! The maxvalue is determined as (2**bitpix - 1).
  ! The specific characteristic of dark images may require some modification
  ! for real-type intensity of image.
  real, parameter :: minvalue = 0.0

  character(len=32768) :: sciname,flatname,outname,buf,key,value, & !mask,
       line,record

!  character(len=666) :: sciname,flatname,outname,mask,buf
!  character(len=15) :: form
!  integer :: icit = 1
!  integer :: i,j, n, nj, mj, nhead
  integer :: i,j, n, nhead, io
  real :: maxvalue
  real :: xmean, xsig
  real, dimension(:,:), allocatable :: FLAT, SCI, OUT
  real, allocatable,dimension(:) :: fbuf
  character(len=80), dimension(:), allocatable :: head

! FITS variables
  integer :: blocksize,naxis,naxes(DIM),pcount,gcount,bitpix,bpix,bitpix0
  logical :: simple,extend,anyf!,exist

  ! input of top line with init values
!  read(*,*,end=90,iostat=j) flatname,mask,icit,bitpix
!  if( j /= 0 ) stop 'Bad format of input data.'

  bitpix0 = 0
  i = 0

  do
     read(*,'(a)',iostat=io,end=20) record
     if( io /= 0 ) stop 'Failed to read of an input.'
  
     call input_rcomment(record,line)

     select case(input_classify(line))

     case(input_key)

        call input_keys(line,key,value)

        if( key == 'Flatframe' ) then
           read(value,*,iostat=io) flatname
           
           ! open flatframe
           call ftopen(20,flatname,0,blocksize,i)
           call ftghpr(20,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,i)
           if( i /= 0 ) goto 666
           allocate(FLAT(naxes(1),naxes(2)),stat=j)
           if( j /= 0 ) stop 'Not enought memory for flat.'
           call ftg2de(20,1,minvalue,naxes(1),naxes(1),naxes(2),FLAT,anyf,i)
           call ftclos(20,i)  
           if( i /= 0 ) goto 666

           ! median of flat
           allocate(fbuf(naxes(1)*naxes(2)),stat=j)
           if( j /= 0 ) stop 'Not enought memory for medianing.'
           fbuf = pack(FLAT,.true.)
           call rmean(fbuf,xmean,xsig)
           deallocate(fbuf)
           write(*,*) 'Flat-frame statistics: mean=',xmean,' dev=',xsig

           ! normalisation
           if( abs(xmean) > epsilon(xmean) )then
              FLAT = FLAT / xmean
           else
              FLAT = huge(0.0)
           endif

           ! allocate buffers for scientific exposure
           allocate(SCI(naxes(1),naxes(2)),OUT(naxes(1),naxes(2)),stat=j)
           if( j /= 0 ) stop 'Not enought memory for scientific image.'

        end if

        if( key == 'Output BITPIX' ) then
           read(value,*,iostat=io) bitpix0
        endif

        if( key == '' .or. value == '' .or. io /= 0 ) then
           write(*,*) 'This record: ',trim(key),trim(value),' parse failed.'
        endif




  ! prepare output mask
!  mj = index(mask,"?")
!  if ( mj > 0 )then
!     do nj = 0, len_trim(mask)
!        if( mask(mj+nj:mj+nj) /= '?' ) exit
!     enddo
!     write(form,"(A,I1.1,A,I1.1,A)") "(I",nj,".",nj,")"
!  endif
!  outname = mask

  ! allocate buffers for scientific exposure
!  allocate(SCI(naxes(1),naxes(2)),OUT(naxes(1),naxes(2)),stat=j)
!  if( j /= 0 ) stop 'Not enought memory for scientific image.'

     case (input_data)
        
        read(line,*,iostat=io) sciname,outname
        if( io /= 0 ) then
           write(*,*) 'Failed to read of filename:',trim(line)
           stop 'Bad format of the input filename.'
        end if

!  do
     ! read name of scientific exposure to subtract
!     read(*,*,end=90,iostat=j) sciname
!     if( j /= 0 ) then
!        write(*,*) "Input scientific name read error. Skipping."
!        goto 665
!     end if
     write(*,'(5a,f0.3,a)',advance="no") '(',trim(sciname),') / (',trim(flatname),'/',xmean,') = '

     ! generate output filename:
     ! by mask
!     if( mj > 0 ) then
!        write(outname(mj:mj+nj-1),form,iostat=i) icit
!        if( i /= 0 .or. index(outname,'*') > 0 ) &
!             stop 'The output filename mask not correspond to sequence number of image.'
!        icit = icit + 1
!     endif

    ! temporary
!     if( mask == '.' ) then
!        outname = sciname
!     endif

     ! work in current dir
!     if( mask == '$' ) then
!        j = index(sciname,'/',.true.) + 1
!        outname = sciname(j:)
!     end if

     if( .not. allocated(FLAT) ) &
          stop 'Flat frame not loaded yet.'

     ! scientific exposure initialisation
     call ftopen(25,sciname,0,blocksize,i)
     call ftghpr(25,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,i)
     if( i /= 0 ) goto 666
     if( naxes(1) /= size(FLAT,1) .or. naxes(2) /= size(FLAT,2) )then
        write(*,*) "The image and flatframe size does not match. Skipping."
        call ftclos(25,i)
        goto 665
     endif

     call ftghps(25,nhead,n,i)
     allocate(head(nhead),stat=j)
     if( j /= 0 ) stop 'Not enought memory for head.'
     do j = 1, nhead
        call ftgrec(25,j,head(j),i)
     enddo

     call ftg2de(25,1,minvalue,naxes(1),naxes(1),naxes(2),SCI,anyf,i)
     call ftclos(25,i)
     if( i /= 0 ) goto 666

     ! flat correction
     where( FLAT /= 0.0 )
        OUT = SCI / FLAT
     elsewhere
        OUT = huge(1.0)
     end where

!     if( bitpix == 0 ) bitpix = bpix

        if( bitpix0 == 0 ) then
           bitpix = bpix
        else
           bitpix = bitpix0
        end if

     ! range checking
     if( bitpix > 0 )then
        maxvalue = 2**bitpix - 1
        where( OUT < minvalue )
           OUT = minvalue
        endwhere
        where( OUT > maxvalue )
           OUT = maxvalue
        endwhere
     end if

     ! check if output file exist, do not overwrite
!     Inquire(file=outname,exist=exist,iostat=j)
!     if( j /= 0 ) stop 'Status of file check failed.'
!     if( exist )then
        ! remove file to overwrite
!        if( mask == '.' .or. &
!             ( mask(1:4) == 'fout' .and. mask(10:13) == 'fits' ) ) then
!           open(26,file=outname,iostat=j)
!           close(26,status='delete',iostat=j)
!           if( j /= 0 ) then
!              write(*,*)"Error occured during remove of output file. Skipping."
!              goto 665
!           endif
!        else
!           write(*,*) "The output file exist (",trim(outname),"). Skipping."
!           goto 665
!        end if
!     endif

     ! write out corrected image
     call ftinit(26,outname,blocksize,i)
     do j = 1, nhead
        call ftprec(26,head(j),i)
     enddo
 !    if( bitpix /= bpix ) then
        call ftukyj(26,'BITPIX',bitpix,'number of bits per data pixel',i)
!        if( bitpix < 0 ) then
!           call ftdkey(26,'BSCALE',i)
!           call ftdkey(26,'BZERO',i)
!        else
        if( bitpix > 0 ) then
           call ftukyj(26,'BSCALE',1,'',i)
           call ftukyj(26,'BZERO',2**(bitpix-1),'',i)
        endif
!     endif
     write(buf,'(f0.3)',iostat=i) xmean
     if( i /= 0 ) buf = ' - '
     call ftphis(26,'Flatframe corrected: '//trim(flatname)//' /'//trim(buf),i)
     call FTPCOM(26,version,i)
     call ftp2de(26,1,naxes(1),naxes(1),naxes(2),OUT,i)
     call ftclos(26,i)
     if( i /= 0 ) goto 666

     write(*,*) '(',trim(outname),')'

665  continue
     if( allocated(head) ) deallocate(head)

  end select

end do
20 continue

!  deallocate(FLAT,SCI,OUT)

666 continue

  if( allocated(FLAT) ) deallocate(FLAT)
  if( allocated(SCI) ) deallocate(SCI)
  if( allocated(OUT) ) deallocate(OUT)

  if( i /= 0 ) then
     write(*,*) 'The following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif

end program flatbat
